#### installing packages if not installed ####

list.of.packages <- c('tidyverse',
                      'ggplot2',
                      'lubridate',
                      'estimatr',
                      'ggthemes',
                      'readr',
                      'readxl',
                      'rdrobust',
                      'gridExtra')
new.packages <- list.of.packages[!(list.of.packages %in% installed.packages()[,"Package"])]
if(length(new.packages)) install.packages(new.packages)

#### attaching libraries ####

suppressPackageStartupMessages(
  
  {
    
    library(tidyverse)
    library(ggplot2)
    library(lubridate)
    library(estimatr)    
    library(ggthemes)
    library(readr)
    library(readxl)
    library(rdrobust)
    library(gridExtra)
    
  }
  
)

#### load data #### 

# loading 911 call data

load(file = "data/calldata.RData")

call = 
  call %>% 
  mutate(blm = ifelse(date >= as.Date("2020-05-29"), 1, 0)) %>% 
  mutate(initiate_off = ifelse(call_type == "-" | call_type == "ONVIEW" | 
                                 call_type == "SCHEDULED EVENT (RECURRING)" | 
                                 call_type == "ALARM CALL (NOT POLICE ALARM)" | 
                                 call_type == "PROACTIVE (OFFICER INITIATED)" | 
                                 call_type == "POLICE (VARDA) ALARM", 1, 0),
         initiate_civ = 
           ifelse(call_type == "911" | 
                    call_type == "TELEPHONE OTHER, NOT 911" | 
                    call_type == "TEXT MESSAGE" | 
                    call_type == "HISTORY CALL (RETRO)" | 
                    call_type == "IN PERSON COMPLAINT", 1, 0))

# call dataset

call_cct_off = call %>% 
  filter(initiate_off == 1) %>% 
  mutate(count = 1) %>% 
  group_by(date) %>% 
  summarize(count = sum(count),
            blm = mean(blm)) %>% 
  mutate(count_calls = count)

call_cct_off$blmrv = 
  call_cct_off$date - min(call_cct_off$date[call_cct_off$blm == 1])


call_cct_noff = call %>% 
  filter(initiate_off == 0) %>% 
  mutate(count = 1) %>% 
  group_by(date) %>% 
  summarize(count = sum(count),
            blm = mean(blm)) %>% 
  mutate(count_calls = count)

call_cct_noff$blmrv = 
  call_cct_noff$date - min(call_cct_noff$date[call_cct_noff$blm == 1])

# calls: 911, los angeles

load("data/calls_la.RData")

calls_noff_ts = lacallsf %>% 
  mutate(count = 1) %>%
  filter(off_call == 0) %>% 
  group_by(date) %>% 
  summarize(calls_noff = sum(count, na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE)) %>% 
  filter(date >= as.Date("2019-01-01"))

calls_noff_ts$blmrv = 
  calls_noff_ts$date - min(calls_noff_ts$date[calls_noff_ts$blm == 1])


calls_off_ts = lacallsf %>% 
  mutate(count = 1) %>%
  filter(off_call == 1) %>% 
  group_by(date) %>% 
  summarize(calls_off = sum(count, na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE)) %>% 
  filter(date >= as.Date("2019-01-01"))

calls_off_ts$blmrv = 
  calls_off_ts$date - min(calls_off_ts$date[calls_off_ts$blm == 1])

#### analysis #### 

# generating covars 

call_cct_off$call_diff = call_cct_off$count_calls - call_cct_noff$count_calls
call_cct_off$call_ratio = call_cct_off$count_calls / call_cct_noff$count_calls

calls_off_ts$call_diff = calls_off_ts$calls_off - calls_noff_ts$calls_noff
calls_off_ts$call_ratio = calls_off_ts$calls_off / calls_noff_ts$calls_noff

# estimates 

rdmod1 = with(call_cct_off, rdrobust(x = blmrv, y = call_diff, p = 1, kernel = "uni"))
rdmod2 = with(call_cct_off, rdrobust(x = blmrv, y = call_ratio, p = 1, kernel = "uni"))
rdmod3 = with(calls_off_ts, rdrobust(x = blmrv, y = call_diff, p = 1, kernel = "uni"))
rdmod4 = with(calls_off_ts, rdrobust(x = blmrv, y = call_ratio, p = 1, kernel = "uni"))

# plot 

plot_df_out = data.frame(
  
  est = c(rdmod1$coef[1], rdmod2$coef[1], rdmod3$coef[1], rdmod4$coef[1]),
  se = c(rdmod1$se[3], rdmod2$se[3], rdmod3$se[3], rdmod4$se[3]),
  outcome = rep(c("Difference", "Ratio"), 2),
  city = c(rep("Seattle", 2), rep("Los Angeles", 2))
  
)

call_diffp1 = plot_df_out %>% 
  filter(outcome == "Difference") %>% 
  ggplot() + 
  geom_hline(yintercept = 0) + 
  geom_point(aes(x = city, y = est)) + 
  geom_errorbar(aes(x = city, 
                    ymin = est - 1.96 * se,
                    ymax = est + 1.96 * se),
                width = 0,
                size = .4) + 
  geom_errorbar(aes(x = city, 
                    ymin = est - 1.645 * se,
                    ymax = est + 1.645 * se),
                width = 0,
                size = .6) + 
  labs(x = "City", y = "RDiT Coefficient (BLM Protest)",
       title = "A. Officer-Civilian Call Difference") + 
  theme_tufte(base_size = 8) + 
  theme(plot.title = element_text(size = 7))

call_diffp2 = plot_df_out %>% 
  filter(outcome == "Ratio") %>% 
  ggplot() + 
  geom_hline(yintercept = 0) + 
  geom_point(aes(x = city, y = est)) + 
  geom_errorbar(aes(x = city, 
                    ymin = est - 1.96 * se,
                    ymax = est + 1.96 * se),
                width = 0,
                size = .4) + 
  geom_errorbar(aes(x = city, 
                    ymin = est - 1.645 * se,
                    ymax = est + 1.645 * se),
                width = 0,
                size = .6) + 
  labs(x = "City", y = "RDiT Coefficient (BLM Protest)",
       title = "B. Officer/Civilian Call Ratio") + 
  theme_tufte(base_size = 8) + 
  theme(plot.title = element_text(size = 7))

#### plots #### 

# seattle

descp1 = call_cct_off %>% 
  ggplot() + 
  geom_point(aes(x = date, y = count_calls),
             alpha = .05,
             size = .1) + 
  geom_smooth(aes(x = date, y = count_calls, group = blm),
              se = FALSE,
              col = "black",
              size = .35) + 
  geom_vline(xintercept = as.Date("2020-05-29"),
             linetype = 2) + 
  labs(x = "Date", y = "Officer Calls",
       title = "C. Officer Calls (Seattle)") + 
  ylim(0, 950) +  
  theme_tufte(base_size = 8) + 
  theme(plot.title = element_text(size = 7))

descp2 = call_cct_noff %>% 
  ggplot() + 
  geom_point(aes(x = date, y = count_calls),
             alpha = .05,
             size = .1) + 
  geom_smooth(aes(x = date, y = count_calls, group = blm),
              se = FALSE,
              col = "black",
              size = .35) + 
  geom_vline(xintercept = as.Date("2020-05-29"),
             linetype = 2) + 
  labs(x = "Date", y = "Civilian Calls",
       title = "D. Civilian Calls (Seattle)") + 
  ylim(0, 950) +
  theme_tufte(base_size = 8) + 
  theme(plot.title = element_text(size = 7))

# los angeles

descp3 = calls_off_ts %>% 
  ggplot() + 
  geom_point(aes(x = date, y = calls_off),
             alpha = .1,
             size = .1) + 
  geom_smooth(aes(x = date, y = calls_off, group = blm),
              se = FALSE,
              col = "black",
              size = .35) + 
  geom_vline(xintercept = as.Date("2020-05-29"),
             linetype = 2) + 
  labs(x = "Date", title = "A. Officer Calls (Los Angeles)",
       y = 'Officer Calls') + 
  ylim(600, 4000) + 
  theme_tufte(base_size = 8) + 
  theme(plot.title = element_text(size = 7))

descp4 = calls_noff_ts %>% 
  ggplot() + 
  geom_point(aes(x = date, y = calls_noff),
             alpha = .1,
             size = .1) + 
  geom_smooth(aes(x = date, y = calls_noff, group = blm),
              se = FALSE,
              col = "black",
              size = .35) + 
  geom_vline(xintercept = as.Date("2020-05-29"),
             linetype = 2) + 
  labs(x = "Date", title = "B. Civilian Calls (Los Angeles)",
       y = "Civilian Calls") + 
  ylim(600, 4000) + 
  theme_tufte(base_size = 8) + 
  theme(plot.title = element_text(size = 7))

#### long term effects #### 

# setting up loop

call_cct_list = 
  list(call_cct_off, call_cct_noff, calls_off_ts, calls_noff_ts)

call_cct_outcomes = 
  c('count_calls', 'count_calls',
    'calls_off', 'calls_noff')

call_cct_list_out = as.list(rep(NA, length(call_cct_list)))

for (i in 1:length(call_cct_list)) {
  
  print(paste0("I iteration ", i))
  
  bw_matrix = matrix(NA, ncol = 5, nrow = length(1:100))
  
  for (l in 1:length(1:100)) {
    
    print(paste0("L iteration ", l))
    
    dftouse = 
      call_cct_list[[i]] %>% 
      as.data.frame %>% 
      filter(blmrv <= -1 | blmrv >= l)
    
    dftouse$blmrv = 
      ifelse(dftouse$blmrv >= 0, dftouse$blmrv - l, dftouse$blmrv)
    
    out = rdrobust(y = as.numeric(dftouse[, call_cct_outcomes[i]]),
                   x = dftouse$blmrv, 
                   p = 1,
                   kernel = 'uni')
    
    bw_matrix[l, 1] = out$coef[3]
    bw_matrix[l, 2] = out$se[3]
    bw_matrix[l, 3] = out$pv[3]
    bw_matrix[l, 4] = call_cct_outcomes[i]
    bw_matrix[l, 5] = (1:100)[l]
    
  }
  
  call_cct_list_out[[i]] = 
    bw_matrix %>% 
    as.data.frame %>% 
    `colnames<-` (c("est", "se", "pv", "outcomes", "days_cut")) %>% 
    mutate(est = as.numeric(as.character(est)),
           se = as.numeric(as.character(se)),
           pv = as.numeric(as.character(pv)),
           outcomes = as.character(outcomes),
           days_cut = as.numeric(as.character(days_cut)))
  
}

# plots 


ltefxp1 = bind_rows(call_cct_list_out[[1]] %>% mutate(category = "Officer"),
          call_cct_list_out[[2]] %>% mutate(category = "Civilian")) %>% 
  mutate(category = factor(category, levels = c("Officer", "Civilian"))) %>% 
  ggplot() + 
  geom_hline(yintercept = 0) + 
  geom_point(aes(x = days_cut, y = est,
                 col = category),
             position = position_dodge(.2),
             size = .15) + 
  geom_errorbar(aes(x = days_cut, 
                    ymin = est - 1.96 * se,
                    ymax = est + 1.96 * se,
                    col = category),
                position = position_dodge(.2),
                width = 0,
                size = .01) +
  geom_errorbar(aes(x = days_cut, 
                    ymin = est - 1.645 * se,
                    ymax = est + 1.645 * se,
                    col = category),
                position = position_dodge(.2),
                width = 0,
                size = .1) +
  labs(x = "Days Cut", y = "RDiT Coefficient (BLM Protest)",
       title = "B. Officer and Civilian Calls (Seattle)",
       col = "Call Type") + 
  scale_color_grey(start = 0, end = .6) + 
  ylim(-375, 75) + 
  theme_tufte(base_size = 8) + 
  theme(plot.title = element_text(size = 7))

abs((call_cct_list_out[[1]]$est - call_cct_list_out[[2]]$est) /
      sqrt((call_cct_list_out[[1]]$se^2 + call_cct_list_out[[2]]$se^2))) > 1.96

ltefxp2 = bind_rows(call_cct_list_out[[3]] %>% mutate(category = "Officer"),
          call_cct_list_out[[4]] %>% mutate(category = "Civilian")) %>% 
  mutate(category = factor(category, levels = c("Officer", "Civilian"))) %>% 
  ggplot() + 
  geom_hline(yintercept = 0) + 
  geom_point(aes(x = days_cut, y = est, col = category),
             size = .15) + 
  geom_errorbar(aes(x = days_cut, 
                    ymin = est - 1.96 * se,
                    ymax = est + 1.96 * se, col = category),
                position = position_dodge(.2),
                width = 0,
                size = .05) +
  geom_errorbar(aes(x = days_cut, 
                    ymin = est - 1.645 * se,
                    ymax = est + 1.645 * se, col = category),
                position = position_dodge(.2),
                width = 0,
                size = .1) +
  labs(x = "Days Cut", y = "RDiT Coefficient (BLM Protest)",
       title = "A. Officer and Civilian Calls (Los Angeles)",
       col = "Call Type") + 
  scale_color_grey(start = 0, end = .6) + 
  ylim(-1800, 0) + 
  theme_tufte(base_size = 8) + 
  theme(plot.title = element_text(size = 7))

abs((call_cct_list_out[[3]]$est - call_cct_list_out[[4]]$est) /
  sqrt((call_cct_list_out[[3]]$se^2 + call_cct_list_out[[4]]$se^2))) > 1.96

grob_calls1 = 
  arrangeGrob(descp3, descp4, descp1, descp2, ncol = 4)

ggsave(plot = grob_calls1, filename = "pics/bigcall1.png", width = 8, height = 2)

grob_calls2 = 
  arrangeGrob(call_diffp1, call_diffp2, ncol = 2)

ggsave(plot = grob_calls2, filename = "pics/bigcall2.png", width = 4, height = 2)

grob_calls3 = 
  arrangeGrob(ltefxp2, ltefxp1, ncol = 2)

ggsave(plot = grob_calls3, filename = "pics/bigcall3.png", width = 8, height = 2)

#### MAKE REGRESSION TABLE FOR PANELS I (difference) & J (ratio) FOR SI ####

library(xtable)

la_calldiff <- as.data.frame(cbind(rdmod1$coef[1],rdmod1$se[3],rdmod1$pv[3],
                                   rdmod1$N[1],(rdmod1$bws[1]*2),rdmod1$bws[1]))

sea_calldiff <- as.data.frame(cbind(rdmod2$coef[1],rdmod2$se[3],rdmod2$pv[3],
                                    rdmod2$N[1],(rdmod2$bws[1]*2),rdmod2$bws[1]))

la_callrat <- as.data.frame(cbind(rdmod3$coef[1],rdmod3$se[3],rdmod3$pv[3],
                                  rdmod3$N[1],(rdmod3$bws[1]*2),rdmod3$bws[1]))

sea_callrat <- as.data.frame(cbind(rdmod4$coef[1],rdmod4$se[3],rdmod4$pv[3],
                                   rdmod4$N[1],(rdmod4$bws[1]*2),rdmod4$bws[1]))

tab <- as.data.frame(rbind(la_calldiff,sea_calldiff,la_callrat,
                           sea_callrat))

tab <- tab %>% mutate(V4 = as.character(V4))
tab <- tab %>% mutate_if(is.numeric,round, digits = 2)
tab <- tab %>% mutate(V3 = "0.000")

cities <- as.data.frame(rbind("LA","Seattle","LA","Seattle"))
outcome <- as.data.frame(rbind("Call Difference","Call Difference","Call Ratio","Call Ratio"))

tab <- as.data.frame(cbind(cities,outcome,tab))
colnames(tab) <- c("City","Officer:Civilian Calls","Coeff","SE","P-Val","N-Val","Effective N","Bandwidth Est.")

xt <- xtable(
  x = tab,
  align = "lllcccccc",
  type = "html",
  caption = "RDiT Coefficients Characterizing Changes in Officer & Civilian Initiated Calls Following the BLM Protest.",
  label = "tab:offVciv_calls_reg_tab"
)


addtorow <- list()
addtorow$pos <- list(nrow(xt))
addtorow$command <- as.vector(c(paste("\\hline\n\\multicolumn{8}{l}{\\emph{All estimates are specified with a uniform kernel and polynomial degree equal to 1.}}\\\\\n",
                                      paste("\\multicolumn{8}{l}{\\emph{Standard errors are robust.}}\\\\\n",
                                            sep=""))))



print(xt,
      format.args = list(big.mark = ","),
      caption.placement = "top",
      include.rownames = FALSE,
      add.to.row = addtorow,
      file = "fig3_IJ_reg_tab.tex",
      hline.after = c(-1, 0, 4)
)





